use strict;
use warnings;
use Data::Dumper;
use Benchmark qw( cmpthese ) ;

my @masas_aas;
use Core::IntConsts qw(get_masas_aa get_masas_mod);
my %masas_aa = %{&get_masas_aa};
my %masas_mod = %{&get_masas_mod};



#Ascore 548
sub prever_fragmentos{
	#Predict the spectrum of a sequence. Allows phosphorilations on
	# STY (+80:s,t,y), water loss on ST (-18:u,v), and Met oxidation (+18:m).
	#It calculates the y and b series for the charge state,
	#ARGUMENTS:
	#     - $_[0]-> sequence
	#     - $_[1]-> charge state (1,2 or 3)
	#     - $_[2]-> mz upper limit for fragment generation
	#RETURN: a hash with {label}->mass.
	my $peptide = $_[0];
	my $z = $_[1];
	my $limit_mz_sup = $_[2];
	my $limit_mz_inf = $_[3];
	my (@aas, @mods);
	my @masas_aas;

	while($peptide =~ /\w/){
		my $aa;
		my $mod;
		if($peptide =~ s/^(\w)\((\w+\,\w+)\)//){
			$aa = $1;
			$mod = $2;
		}
		elsif($peptide =~ s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
		}
		else{
			$peptide =~ s/^(\w)//;
			$aa = $1;
			$mod = '';
		}
		push @aas, $aa;
		push @mods, $mod;
	}

    for(0..$#aas){
		if($mods[$_]){
			if($mods[$_] =~ /\,/){
				my ($m1, $m2) = split /\,/, $mods[$_];
				$masas_aas[$_] = $masas_mod{$m1} + $masas_mod{$m2} + $masas_aa{$aas[$_]};
			}
			else{
				$masas_aas[$_] = $masas_mod{$mods[$_]} + $masas_aa{$aas[$_]};
			}
		}
		else{
			$masas_aas[$_] = $masas_aa{$aas[$_]};
		}
	}
#sigue..................
}


#Ascore 318
sub peptide_collection{
	my $pept = $_[0];
	#Limitado a 1/dos chemtags por aa o 1 chemtag + 1 phos
	#			2/Los chemtag no cambian de sitio
	my @posibles_fosfo_sites;
	my @posibles_deshidrat_sites;
	my $num_fosfos = 0;
	my $num_deshidrats = 0;
	my $num_chemtag = 0;
	my $alrdy_one;

	my (@aas, @mods);
	my $aa_num = 0;
	while($pept){ #=~ /\w/
		my ($aa, $mod);
		if($pept =~ s/^(\w)\((\w+),(\w+)\)//){
			$alrdy_one = 0;
			$aa = $1;
			my @mod1 = ($2, $3);
			push @aas, $aa;
			#Calculo el numero de fosfos/chemtags/DH
            foreach(@mod1){
                if($_ == 21){
                    $num_fosfos++;
                }
                elsif($_ == 23){
                    $num_deshidrats++;
                }
                else {
					#Limitado a dos chemtags por aa
                    if ($alrdy_one){
						my $frst = pop @mods;
						push @mods, "$frst,$_";
					} else {
						push @mods, $_;
						$alrdy_one = 1;
					}
                    if($aa_num > 0){        #Si hay chemtag en un aa no N-terminal,
                        $aa_num++;          #no se puede modificar con nada mas
                        next;               #ojo este next tiene que ser un label
                    }
                }
            }
		}
        elsif ($pept =~ s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
			#Calculo el numero de fosfos/chemtags/DH
			if ($mod == 21){
				$num_fosfos++;
			}
			elsif ($mod == 23){
				$num_deshidrats++;
			}
			elsif ($mod == 214 or $mod == 737){
				$num_chemtag++;
				if ($aa_num > 0){        #Si hay chemtag en un aa no N-terminal,
					$aa_num++;           #no se puede modificar con nada mas
					next;
				}
			}
        }
		else{
			$pept =~ s/^(\w)//;
			$aa = $1;
			$mod = '';
			push @aas, $aa;
			push @mods, $mod;
		}

		if($aa =~ /[S,T,Y]/){
			push @posibles_fosfo_sites, $aa_num;
		}
		if($aa =~ /[S,T]/){
			push @posibles_deshidrat_sites, $aa_num;
		}
		$aa_num++;
	}
	#sigue..........
}


my $pep = "L(214)QDSSDPDTGSEEEGS(23)SRLSPPHSPR";
my $z = 1;

#Integration 375
sub calculate_mz{
	#my $pep = $_[0];
	#my $z = $_[1];
	my $nt;
	my @aas;
	my @mods;

    while($pep =~ /\w/){
		if($pep =~ s/^(\w)\((\w+),(\w+)\)//){
            push @aas, $1;
			push @aas, 'Not';
            push @mods, $2;
			push @mods, $3;
		}
		elsif($pep =~ s/^(\w)\((\w+)\)//){
            push @aas, $1;
            push @mods, $2;
		}
		else{
			$pep =~ s/^(\w)//;
			push @aas, $1;
			push @mods, 'Not';
		}
	}
	my $mass = 18.010565;      #H2O
	for (0..$#aas){
		$mass = $mass + $masas_mod{$mods[$_]} + $masas_aa{$aas[$_]};
	}

	return ($mass + $z * 1.007825) / $z;
}


#New 300% faster, Integration 375
sub calculate_mz_2{
	#my $pep = $_[0];
	#my $z = $_[1];
	my $nt;
	my $mass = 18.010565;      #H2O

    while($pep =~ /\w/){
		if($pep =~ s/^(\w)\((\w+),(\w+)\)//){
            $mass +=   $masas_aa{$1} + $masas_mod{$2} + $masas_mod{$3};
		}
		elsif($pep =~ s/^(\w)\((\w+)\)//){
            $mass +=   $masas_aa{$1} + $masas_mod{$2};
		}
		else{
			$pep =~ s/^(\w)//;
			$mass +=   $masas_aa{$1};
		}
	}

	return ($mass + $z * 1.007825) / $z;
}

#"L(214)QDSSDPDTGSEEEGS(23)SRLSPPHSPR" => '2892.32435',
#print Dumper calculate_mz('L(214)QDSSDPDTGSEEEGS(23)SRLSPPHSPR', 1);
#print Dumper calculate_mz_2('L(214)QDSSDPDTGSEEEGS(23)SRLSPPHSPR', 1);
#
#__END__

cmpthese(500000, {
        a => \&calculate_mz,
        b => \&calculate_mz_2,
        #c => \&pep,
        #d => \&pep
        } );



1;
